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Abstract 

We consider an ensemble of globally coupled phase oscillators whose interac- 
tion is transmitted at finite speed. This introduces time delays, which make 
the spatial coordinates relevant in spite of the infinite range of the interac- 
tion. We show that one-dimensional arrays synchronize in an asymptotic state 
where all the oscillators have the same frequency, whereas their phases are 
distributed in spatial structures that -in the case of periodic boundaries- can 

propagate, much as in coupled systems with local interactions. 
PACS: 05.45.Xt, 05.65.+b 
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I. INTRODUCTION 



Standard models for studying collective complex behavior in natural systems consist 
typically of ensembles of interacting dynamical elements Jl| . Such kind of models has proven 
to be extremely versatile in the mathematical description, both analytical and numerical, 
of a wide variety of phenomena within the scopes of physics, biology, and other branches of 
science 0. According to the range of the involved interactions, these models can be divided 
into two well distinct classes. Local interactions -which are paradigmatically represented in 
reaction-diffusion systems [T]]- give rise to macroscopic evolution in which space variables 
play a relevant role, such as spatial structures and propagation phenomena. On the other 
hand, with global interactions -where the coupling range is of the order of, or larger than, 
the system size- space becomes irrelevant and collective behavior is observed to develop in 
time, typically, in the form of synchronization ||. 

An essential model of globally coupled elements is given by a set of N identical oscillators 
described, in the so-called phase approximation, by phase variables <fii(t) {i = 1,...,N). 
Their evolution is governed by the equations 

e N 

(p i = uj + —J2 sin ( ( Pj- ( Pi)- (!) 

iV i=i 

It is known that, for any value of the coupling intensity e, all the elements converge to a 
single orbit whose frequency u coincides with that of an individual oscillator ||. In this 
case, e~ l measures the time required to reach such synchronized state. 

In this note we present results on a generalization of the above model, where time delays 
are introduced. The effect of time delays in synchronization phenomena has already been 
considered for two-oscillator systems, both periodic |4| and chaotic ||. Ensembles with lo- 
cal interactions and globally interacting inhomogeneous systems have also been studied 
0. None of these contributions make however explicit reference to the relevant case where 
interactions are global but their propagation occurs at a finite velocity v. This situation, 
which naturally introduces time delays, provides a realistic description of highly connected 
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systems where the time scales associated with individual evolution and with mutual signal 
transmission are comparable. Instances of such systems are neural and informatic networks 
||, and biological populations with relatively slow communication media -such as sound 
||. Our main result is that, since a finite signal velocity makes spatial variables relevant 
even when interactions are global, globally coupled ensembles with time delays exhibit typ- 
ical features of systems driven by local interactions, in particular, structure formation and 
propagation. 

II. THE MODEL AND ITS SOLUTION FOR SHORT DELAYS 

We consider an ensemble of N identical oscillators in the phase approximation, governed 
by the equations 

e N 

Ut) = u + - J2 M<f>i(t - Tij) - Ut)l (2) 

iV 3=1 

where Ty = d^jv is the time required for the signal to travel from element j to element i at 
velocity v, and d^ is the distance between i and j. Note that coupling is still global, since 
its intensity e does not depend on the distance between elements. However, the relative 
position of the oscillators becomes now relevant through time delays. 

The full specification of our system requires to fix the topology and the metric properties 
of the ensemble, by fixing the values dij for all i,j — l,...,N. Moreover, initial conditions 
for <pi must be provided. In the case of delay equations like (|J), it is necessary to specify 
the evolution of 0j at times prior to t = up to a time Tj = — max{rjj}j [|IIJ]. In the 
following we shall assume that for t < the oscillators evolve independently from each other 
at their proper frequency u and with random relative phases. Namely, for t < we have 
<f>i(t) = ut + 0j(O), where 0j(O) is chosen at random from a uniform distribution in [— n, ir). 
At t — coupling in switched on, so that we formally have a time-dependent coupling 
intensity e(t) = e9(t), where 6 is the Heaviside step function. 

Through extensive numerical calculations for a variety of topologies, ranging from one- 
dimensional arrays to tree (ultrametric) structures, we have found that the system evolves 
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to a state where all the oscillators have the same frequency. On the other hand, in contrast 
with the case without time delays ||, their phases can be different. This asymptotic state 
corresponds thus to a situation of frequency synchronization. The long-time evolution of 
each oscillator can then be written as <pi(t) = Qt + ipi, where in general ipi 7^ ipj for i 7^ j. 
The fact that these phases are different could have been expected for topologies where not 
all the elements are equivalent -for instance, when boundaries are present. As we show 
later, however, such states are also found in homogeneous topologies. In this case, they are 
associated with propagating structures. 

In general, the synchronization frequency is different from the proper frequency of each 
oscillator, Q 7^ uj. According to (fj), the synchronization frequency satisfies 

e N 

Q = to - — sin(fir ij - ipj + ipi). (3) 

Note that the sums Si = J2j sin(fl7ij — ipj + ipi) are in general different for each i. However, 
their numerical value must coincide if the synchronization frequency is to be well defined. 
For a given value of Q, this constraint provides N — 1 independent equations for the phases 
ipi- 

Sl = s 2 = ... = S N . (4) 

Since phases are defined up to an additive constant we can choose for instance ip\ = 0, and 
solve these equations for ip2, ■ ■ - ipN- Then, Q can be found self-consistently from (|^). For 
large values of N, and due to the involved nonlinearities, this results to be a quite hard 
numerical problem. 

An approximate solution can however been found in the case of short delays, i.e. close to 
the situation where the system is also synchronized in phase, ipi = ipj for all Assuming 
that \ flTij — ip j + ipi I 1 , we can write 

s, « xxfirtf - iPj + iPi) = Nn( n ) -$>,• + Nip h (5) 
3 3 

where (rj) = A^ -1 J2j T %j is the average of the time delays associated with element i. Taking 
into account Eq. (f|), we get 
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ViWtf-ftfo). (6) 

where \I/ is a constant, independent of z, that can be chosen arbitrarily. This result indicates 
that, in this short-delay limit, oscillators with small average delays are relatively ahead in 
the evolution, as their phases are larger. This is plausibly due to the fact that, in average, 
they receive the information on the system state before other elements with larger values 
of (Tj), which are thus relatively retarded. Note moreover that in a homogeneous topology 
all the elements are equivalent, so that (Tj) is the same for all oscillators. In this case, the 
system is also synchronized in phase. 

Within the approximation of short delays, the synchronization frequency is given by 

where ((r)) = A r_1 J2i( T i) is the overall average delay. It therefore results that Q is smaller 
than the proper frequency of each oscillator. 



III. ONE-DIMENSIONAL ARRAYS 

In this note, we specifically focus the attention on the case of one-dimensional arrays. 
Two different topologies are considered, namely, with periodic boundary conditions -where 
all the elements are equivalent- and with free boundaries -where the neighborhood of each 
element depends on its distance to the center of the array. For periodic boundary conditions, 
which we consider first, the distance between two elements is not univoquely defined, since 
it can be measured around the ring in both directions. We fix by taking the minimum 
of these values, namely, c?y = min{|z — j\, N — \i — The delay time is thus = 
Tomin{|i — j\,N — \i — j\}, where tq is the time required for the signal to travel between 
nearest neighbors. 

In equations (0), the proper frequency uj can be used to define time units so that, without 
loss of generality, we fix u = 1. Moreover, our numerical simulations are restricted to the 
case e = 1. As a matter of fact, we have found that other coupling intensities do not produce 
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qualitatively different results. Note that this would not be the case if the oscillators had 
chaotic individual dynamics. In such situation, the value of e is expected to control the 
existence of synchronized states. 

We have solved numerically equations (0) for ensembles of N = 10 2 to 10 4 oscillators 
with a standard finite-difference scheme. For small values of r we find that the above 
described random-phase initial conditions evolve to a state of synchronized frequency where 
the phases of all oscillators coincide, ipi = ipj for all i, j. This fully synchronized state is 
completely analogous to that of globally coupled identical oscillators without time delays, 
and corresponds to the approximate solution @ for the present homogeneous topology. 
In this case, (|3]) becomes an autonomous equation for Q. The sum in the right-hand side 
can in fact be explicitly evaluated -though its expression depends on N being even or 
odd- and the synchronization frequency can be found numerically by standard methods. In 
general, this equation admits several solutions. For the values of tq where the state of phase 
synchronization is encountered, however, there is only one possible value of Q. 

At t ~ 5A^ _1 a qualitative change occurs. Above this critical value, the asymptotic 
synchronized state is not characterized by a homogeneous phase anymore. Instead, the phase 
varies linearly along the system, in such a way that a phase difference Aip = ±2tt accumulates 
in a whole turn around. The sign of Aip is defined by the initial condition. Symmetry 
considerations, in fact, indicate that both signs will be found with equal probability over 
the set of initial conditions that lead to this kind of asymptotic state. The individual phases 
are given by ipi = ipo + iSip, with 5ip = ±2tt/N and ipo an arbitrary constant. Due to the 
time evolution of (fii(t) = Qt + i5ip + 4>o a structure propagates around the system at velocity 

vi = -n/5ip. 

Similar qualitative changes are found at larger values of r . For r ~ HAT -1 , 16A^ _1 , . . ., 
the asymptotic states modify their phase structure in such a way that the phase difference 
around the whole system, mAip = ±2irm with m = 2,3, . . ., increases progressively. The 
corresponding individual evolution is <pi(t) = Qt + imSip + "ipo, which defines a propagation 
velocity V m = —Q/m5ip. The synchronization frequency is given by 
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Q = uj — — sin[fir min{|i — j\,N — \i — j\} + (i — j)m5ip}- (8) 
i 

For m = this reduces to the case of full synchronization found for small To- Figure [I] 
shows the solutions of equation (|8|) for various values of m, and N = 100. Bolder curves 
indicate the intervals where each mode has been observed in the numerical calculations with 
random-phase initial conditions. Note the zones where more than one solution exist for 
m = and m = 1. 

Are the transitions observed at the above quoted values of r actual bifurcations, asso- 
ciated with changes in the stability of the asymptotic states? In view of the difficulty of 
dealing with the linear stability problem for a many-dimensional system with time delays 
such as (HD [pLOU , we choose to answer this question by numerical means. For a given value 
of To we calculate the frequency Q of a given mode m from equation (||D and generate an 
initial condition which corresponds to that mode added with a certain -typically random- 
small perturbation. Then, we run the evolution and study the asymptotic behavior. This 
has been carried out for m — 0, . . . , 3 at several values of tq in (0, 0.2), for a 100-oscillator 
ensemble. In almost all cases, it has been found that for sufficiently small perturbations the 
considered states are stable for any value of To. The only exceptions seem to be the states 
whose frequencies are multiple solutions of equation since in this case the only stable 
state correspond to the smallest frequency. 

The observed transitions are therefore not related to stability changes in the propagation 
modes. Rather, several modes coexist and the system is multistable. The specific asymptotic 
state is thus selected by the initial condition. The fact that from the random-phase initial 
conditions considered previously the system evolves to a well defined synchronous mode, 
whose order m grows with r , suggests that the attraction basins of the various solutions 
could considerably vary in size as To changes. Indeed, from a probabilistic viewpoint, most 
initial conditions are of the random-phase type. Initial conditions that, for a given value 
of Tq, do not evolve to the mode marked with a bold line in Fig. p] should be considered 
probabilistically rare. 
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We consider now the case of a one-dimensional array with free boundaries. Here, the 
distance between elements can be defined in the standard form, djj = \i — j\, so that the 
time delays are r%j = tq\i — j\. In this topology sites are not equivalent. Delays for elements 
near the center of the array are in average lower than for elements towards the boundaries. 
As a consequence, no homogeneous stable states are expected for the coupled ensemble. 
Numerical results show that, in fact, in the asymptotic evolution all the oscillators have the 
same frequency, given by 

= w - ^^sin(fir |i -j\ - tjjj + tpi), (9) 

3 

but ipi ipj if i 7^ j for any nearest-neighbor time delay tq. Unexpectedly, however, the 
associated spatial structures not always preserve the topological symmetry of the system, 
as shown in the following. 

Our numerical calculations for the case of free boundaries correspond to a 100-oscillator 
ensemble with the random-phase initial conditions described above. For small values of r 
we find a symmetric asymptotic pattern, ipi = ipN/2-i, where the central elements have larger 
phases than near the boundaries (Fig. |2| for r = 0.02). This structure corresponds to the 
approximate solution (|6|) which, for this topology, predicts a parabolic phase profile with 
a maximum at the center of the array. Beyond a critical value tq m 0.025 random-phase 
initial conditions are instead attracted towards an asymmetric structure, where the phase 
varies in — ~ ^ from one end to the other, and attains a maximum in between. 
Figure shows such structure for tq = 0.05. In average, of course, half of the realizations 
produce the symmetric counterpart of this asymptotic state. The situation changes again at 
To ~ 0.06. Beyond this point, stationary structures are again symmetric, as shown in Fig. 
for t = 0.1. They result however to be more complicated than for small r , with inflection 
points at i w N/4 and a much flatter maximum. A new critical point occurs at r ~ 0.11, 
where phase structures become asymmetric once more (see Fig. |2], for tq = 0.12). The phase 
variation between the ends is similar to that observed for smaller r but the intermediate 
geometry is considerably more complex. 
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An analytical or semi-analytical study of these structures -including their existence and 
stability properties- requires considering the consistency problem discussed in connection 
with Eq. (|3p. Fixing ipx = 0, the N — 1 equations for ipi (i = 2, . . . , N) read here 

^sin(fir |i - j\ - ipj +ipi) = ^sin[Or (j - 1) - 4>j}. (10) 
j j 

This problem will be discussed in detail in a separate publication [|lT |. Let us stress for 



the moment that, though the appearance of spatial structures was to be expected in an 
inhomogeneous system as the present one- dimensional array with free boundaries, these 
patterns are found to exhibit an unexpected richness upon variation of tq -including, in 
particular, symmetry breaking. 



IV. SUMMARY AND DISCUSSION 

We have found that an ensemble of identical globally coupled oscillators with finite in- 
teraction velocity, which gives origin to time delays, evolves to an asymptotic state where all 
the oscillators have the same frequency but different phases. Generally, the synchronization 
frequency differs from the proper frequency of individual oscillators, so that the dynamics 
of each element in the collective asymptotic motion does not coincide with its individual 
(uncoupled) dynamics. Phases, in turn, are distributed according to spatial patterns with 
nontrivial topological and dynamical properties. Specifically, in a one-dimensional periodic 
array several asymptotic states coexist, corresponding to propagation modes with different 
velocities. In a bounded one-dimensional array we have observed stationary phase structures 
whose symmetry properties depend on the size of time delays. These features, which are 
reminiscent of the behavior of reaction-diffusion systems with local interactions, point out 
sharp differences with the collective motion of coupled oscillators without time delays. 

It is natural to ask whether any structure similar to those described above is observed in 
other, more complex topologies. To advance an answer to this question, we have performed 
a preliminary analysis of a two-dimensional array of 20 x 20 elements with periodic boundary 
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conditions. In this case, each element can be labeled by two indices, i x and i y , according 
to its Cartesian coordinates in the lattice. For algorithmic convenience we have defined the 
distance between elements as dij = = mm{\i x — j x \, L — \i x — j x \} + min{\i y — j y \, L — 

\i y — j y \}, with L = 20 in our case. As above, the delay time is = r dij. In complete 
agreement with the one-dimensional analog, we have here found that for small To the system 
synchronizes both in frequency and phase. Beyond a critical value tq « 0.025, instead, 
propagating phase patterns are observed. Figure |] shows the simplest of these patterns, 
corresponding to the propagation mode with m x = m y = 1. 

Work in progress is being devoted to the detailed characterization of the phase structures 
described in this note, as well as those that could arise in other topologies. The next step 
will be to study the effects of the present kind of time delays in ensembles formed by chaotic 
oscillators, where coupling competes as a stabilizing mechanism against the inherently un- 
stable dynamics of individual elements. 
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FIGURES 
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FIG. 1. Synchronization frequency Q of the asymptotic modes m = 0, ...,3 in a 
one-dimensional ensemble of N = 100 globally coupled oscillators with periodic boundary con- 
ditions, as a function of the nearest-neighbor delay time tq. 
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FIG. 2. Stationary phase patterns in a one-dimensional 100-oscillator array with free bound- 
aries, for various values of To- Without loosing generality, we have fixed V'l = 0. 

12 




FIG. 3. Snapshot of a propagating structure in a two-dimensional 20 x 20-oscillator array with 
periodic boundary conditions. Dark and light zones correspond to phases near zero and ±7T, 
respectively. 
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